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The ground state of the square lattice bilayer quantum antiferromagnet with nearest and 
next-nearest neighbour intralayer interaction is studied by means of the modified spin wave 
method. For weak interlayer coupling, the ground state is found to be always magnetically 
ordered while the quantum disordered phase appear for large enough interlayer coupling. 
The properties of the disordered phase vary according to the strength of the frustration. 
In the regime of weak frustration, the disordered ground state is an almost uncorrelated 
assembly of interlayer dimers, while in the strongly frustrated regime the quantum spin 
liquid phase which has considerable Neel type short range order appears. The behavior of 
the sublattice magnetization and spin-spin correlation length in each phase is discussed. 
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I. INTRODUCTION 



The spin-1/2 square lattice Heisenberg model is now widely believed to have an antiferro- 
magnetic long range order in the ground state. However, it is expected that the strong 
quantum fluctuation in this system may lead to the destruction of the long range order 
with the help of some additional mechanism. In this context, the square lattice antiferro- 
magnetic Heisenberg model with nearest and next-nearest exchange interaction (hereafter 
called Ji — J2 model) P-p!^ and the bilayer Heisenberg model |TP-^ have been studied 
extensively. Both of these models are expected to have the quantum disordered ground 
state for appropriate parameter regime. However, because the mechanisms leading to the 
quantum disordered phase in these two models are of very different nature, it must be most 
interesting to study their interplay in the bilayer Ji — J2 model. 

In the Ji — J2 model, the competition between the nearest neighbour interaction Ji 
and the nearest neighbour interaction J2 introduces the frustration in spin configuration 
which enhances the quantum fluctuation. However, the conclusion about the presence of 
the quantum disordered state in this model is still controvertial even in the most frustrated 
regime. 

On the other hand, in the bilayer model, if the interlayer antiferromagnetic coupling 
is strong enough, the spins on both layers form interlayer singlet pairs and the quantum 
fluctuation is enhanced leading to the quantum disordered state. In other words, the anti- 
ferromagnetic interlayer coupling reduces the effective spin magnitude. Therefore this model 
may be regarded as the physical realization of the single layer Heisenberg model with spin 
less than 1/2. Actually, in the study of the single layer Ji — J2 model, there are consider- 
able number of works which cast doubt on the presence of quantum disordered phase even 
for S' = 1/2 and J2/ Ji = 0.5. f7|-|T0| But some of these works also predict the presence of 
quantum disordered phase for S* < 1/2 which is unreachable within the single layer model. 
I^J^ The bilayer Ji — J2 model can effectively realize such situation. 

This paper is organized as follows: The bilayer Ji — J2 model and its clasical ground 
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state are explained in the next section. In section 3, the modified spin wave approximation 
is apphed to this model. The phase diagram and the behavior of physical quantities are 
presented in section 4. The last section is devoted to summary and discussion. 



II. BILAYER Ji - J2 MODEL 

The Hamiltonian of the bilayer Ji — J2 model is given as follows, 

H = J, {SfSf + SfSf) + J, Y: iStSf + SfSf) + J,j:stSf, (2.1) 

where S!^ is the spin operator with magnitude 5* on the z-th site of the layer a(yU = A 01 B). 
The expression ^ and ^ denote the summation over the intralayer nearest neigh- 

nn nnn 

bour pairs and next nearest neighbour pairs, respectively. The last term represents the 
interlayer coupling. All exchange couplings are assumed to be antiferromagnetic. In the 
following, we denote the ratios J2I J\ = a and J^/ Ji = {3 and take the energy unit Ji = 1. 

In the classical limit, the ground state is the Neel state and the coUinear state according 
as a < 0.5 or a > 0.5. Actually, if the quantum fiuctuation is completely neglected, 
infinite number of ground state configurations are degenerate for a > 0.5. |6|,|25|| However, 
for the single layer Ji — J2 model, it is known that this degeneracy is lifted by the quantum 



fiuctuation and the coUinear ground state is chosen. pSf It is straightforward to extend this 
argument to the bilayer model. Therefore, in the following, we only consider the coUinear 
order for a > 0.5. 

Thus we consider the following two types of spin configuration in the classical limit; 

' : Neel (N-) configuration, (2.2) 

Sf' = 5(_i)™+"+i 



Sf' = ^(-1)"^ 

: coUinear (C-)configuration, (2.3) 

Sf' = S(-l)™+i 
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where the position of the i-th site is denoted by {m,n). In the following, we treat 
the quantum fluctuations around these configurations by means of the modified spin wave 
method. |E3 



III. MODIFIED SPIN WAVE APPROXIMATION 



Based on the classical configurations explained in the last section, let us introduce the 
Dyson-Maleev transformation ||26| , |27[| for each spin as follows: 
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for m + n = odd. 



(3.4) 
(3.5) 

(3.6) 

(3.7) 
(3.8) 



for the N-configuration (|2.2|) , and 
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Sf^ = - 4a.), Sf^ = ^(1 - blk)h 



for m = even, 



> for 171 = odd. 
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(3.9) 
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(3.12) 
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for the C-configuration ( p.3|) . 

For the N-configuration ( p.2|) , the Hamiltonian ( |2.1| ) is rewritten as. 



H = ^ {—25'^ + S{alai + a^aj — aloj — aittj 

<i,j> nn 

+ b\h + b]bj - bib] - b,bj) + a\{a] - a,faj/2 + b\{b] - bi)%/2} 

+ X! {25*^ — S{al — aj)(aj — aj) — a|aj(aj — aj)'^6ie/2 — (a| — a])'^aiaj6io/2 

nnn 

- S{b\ - b])ik - b,) - b\b]{h - b,f5,j2 - {b\ - b])%b,5,j2} 

+ + + ^l^i - 4bl - aA) + alibi - ai)\/2}, (3.14) 

i 

where Sie{Sio) = 1 or according as m + n = even (odd) or odd(even) where Vi = {m,n). 



For the C-configuration (|2.3| ), we have 



H = ^ {—25'^ + S{alai + a^aj — a|aj — ajOj 

nnx 

+ b% + b]b, - bib] - bib,) + al{a] - aifaj/2 + bl{b] - bifbj/2} 



+ ^ {25*^ — S{al — aj)(aj — a^) — ala]{ai — ajY5iey/2 — (a- — a]Yaiaj5ioy/2 

<i,j> nny 

- S{bl - b])ibi - b,) - [bl - b]fhb,5,,y/2 - blb]{b, - b,)H,j2} 



<i,i> 

+ foffoi + b]bj - bib] - bibj) + 4(4 - a,)2a,/2 + bl{b] - bifbj/2} 

+ /^El--^' + -^(^a. + - 4^1 - ci^h) + 4(^1 - ai)'&i/2}, (3.15) 

i 

where 5jey(5io?/) = 1 or according as m = even (odd) or odd(even). The notations E 
and E denote the summation over the intralayer nearest neighbour pairs in x and y 
direction, respectively. 



Following Takahashi, ||2J] we assume the constraint that the sublattice magnetization 



vanish as expected for the two dimensional spin system with continuous symmetry at finite 



temperatures. [p§ 



5 = E < 4«. >= E < > • (3.16) 

i i 

We impose this condition even in the ground state where the long range sublattice mag- 
netization may be present. This means that the average is taken over the direction of the 
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sublattice magnetization even in the ordered phase. Nevertheless, we can calculate the sub- 
lattice magnetization from the long range part of the correlation function which originate 
from the Bose condensate of the bose fields and bi. [^,0 Although the validity of this 



procedure is not well founded, these are the common features of the modified spin wave 
method and we do not discuss this point further. 

We treat the nonlinear terms in ( p. 141) and (|3.15| ) by the mean field approximation as. 



^MF _ ^ {A{alai + a]aj - a\a\ - a^aj + h\hi + h]h, - b\b] - hhj) + 2/\^ - 4SA} 

nn 

nnn 

+ (5Y,{/\AB{a\ai + feffei - affoj - aA) + - 2^Aab} 

i 

- Y.{ii{S - a\a,) + ii{S - h%)}, (3.17) 

i 

for the N-configuration phase. Here yU is the Lagrangian multiplier corrresponding to the 
constraint ( p.lC ). The order parameters are defined by 



< alttj >=< bjbj >= Qxy for rj = ri + d, (3.18) 

< aja] >=< bib] >= A for rj = ri + p, (3.19) 

< alb] >=< a^bi >= Aab, (3.20) 

where p is the vector to the nearest neighbour sites and d to the next nearest sites. 
For the C- configuration, we have 



^MF _ j2 {A{alai + a]aj - a\a] - aiUj + b% + b]bj - bib] - hbj) + - ASA} 
+ 5Z - - «i) - - b]){bi - bj) + ASqy - 2g^} 

+ X! {^xyial^i + — — OjOj + blbi + b]bj — bib] — bibj) + 2A^y — ASA^y} 

nnn 

+ /5^{AAB(alai + b% - albl - aA) + A^b - 2^Aab} 

i 

- T.{^'iS - alai) + ^i{S - blh)]. (3.21) 

i 

Here, the order parameters are defined by 



< alttj >=< h\hj >= Qy for rj = ri + py, (3.22) 

< ala] >=< bib] >= for rj = n + p^, (3.23) 

< a|aj >=< b\b] >= A.^y for Vj = ri + d, (3.24) 

< alb] >=< aibi >= Aab, (3.25) 

where p^ and py is the vector to the nearest neighbour sites in x and y direction, respectively. 
These mean field Hamiltonians are transformed into the fourier space as 

^MF ^ ^ >|^(;^) (^t^^^ ^ ^t_^^_^ ^ ^t^^^ ^ ^t_^^_^) (3 26) 
k 

- r(fc)(a[.a^_^ + afea_fc + + 6fc&_fc) 

+ PAABia^af^ + alf^a_f^ + + blf^b_f^ - a^b^^ - af.b_f. - b^a^f^ - bj.a_j.)} 

+ ^0, (3.27) 

where 

A(fc) = 4(A - aq,y{l - cos{K) cos{ky))) + /i, (3.28) 

r(fc) = 4A7(fc), (3.29) 

7(fc) = ^(cos(A;,.) + cos(A;j,)), (3.30) 

Eo = N{4A^ - 8SA + 4a{4Sq,y - 2qly) (3.31) 

+ /?(AiB - 2^Aab) - 2/i^}, (3.32) 

for the N-configuration and 

A(fc) = -2gj,(l - cos(fcy)) + /i, (3.33) 

r(A;) = 2Aa, cos(A;a;) + 4aA^y cos(A;^) cos(A;y), (3.34) 

= N{{2Al - ASA^ + ASqy - 2gJ + 2a(2A2^ (3.35) 

- 4^A,.j,) + P{Al^ - 2SAab) - 2/i^}, (3.36) 

for the C-configuration. The summation ^ is taken over the left half of the Brillouin zone 

k 

{kx > 0) because the momenta k and —k are explicitly wrtten in ( |3.26| ). The number of the 



lattice sites in each layer is denoted by A^. 



These Hamiltonians are diagonalized as 



k 



where 



E±{k) = J7j{k)^-{r{k)±Sy, S = (3Aab, v{k) = A{k) + 6, 



by the Bogohubov transform 



where 



ch^^ 



\AE±{k) 



+ 1 , sh^^ 



\ 2 [E^{k) 



1 sgn(rfe±5) 



The self consistent equations for the order parameters are 

A 



^AB 



iV^ AE^ik) 

j_ 5±r(fc) ^ 



for the N-configuration and 
A 



'■xy 



A.= 



^AB 



1 ^±r(fc) c 

AT 2^ AC^ IU\ ^ > 



4i?^(fc) 
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for the C-configuration. The summation excludes fc's for which E±{k) = 0. The 

k,± 

quantities Nq and Nq are the amphtudes of the Bose condensate of the bose particles 
represented by the operators af, and Pf. in the Neel phase and collinear phase, respectively. 
They correspond to the magnitude of the Neel or collinear long range order. These 



quantities vanish in the disordered phase. The constraint ( |3.16| ) is rewritten in the form 



^-^^S'^ly + A'o. (3.51) 

in both cases where A^^o stands for A*"^ or A*"^. In order that Nq remains finite, the excitation 
spectrum must have zero modes. This requirement fixes the value of /i in the ordered phase. 
On the other hand, for the solution corresponding to the disordered phase, the energy 
spectrum has no zero mode and the value of /i must be fixed so that the self-consistent 
equations are satisfied with A^o = 0. The ground state energy Eq is given as follows: 

Eg = E^ = Ar(-4A2 + Aaqly - pAl^^) : N-configuration, (3.52) 
Eg = E^ = N{-2Al + 2ql - AaAly - (3A\^) : C-configuration. (3.53) 

If the self consistent equations have more than two solutions, we choose the thermodynam- 
ically stabe solution comparing the ground state energy. 



IV. PHASE DIAGRAM 

We have solved numerically the set of equations for S = 1/2 and found 4 kinds of phases: 

i) Neel ordered phase : A, Aab, <ixy, Nq > 

ii) Collinear ordered phase : A^y, A^., Aab, Qy, > 

iii) Spin hquid phase with Neel-type correlation (NSL phase): A, Aab,Q'xj/ > 0,Nq = 

iv) Interlayer dimer phase (ILD phase): A = q^y = A^y = A^ = qy = Nq = Nq = 
0,Aab = 

The spin fuild phase with collinear-type correlation ( A^y, A^, Aj^-B,<ly > 0,^0 = ) is 
only found as a metastable state. The ground state phase diagram is shown in Fig. |l[ 
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For small /3, there appears no disordered phase. The ground state changes from the Neel 
ordered phase to the collinear ordered phase by the first order transition around 
For larger values of (3, however, there appears a NSL phase around a ~ 0.5. In this phase, 
the short range intralayer singlet order parameter A remains finite, although the long range 
order is absent. The transition between the Neel phase and NSL phase is of the second 
order, while between the collinear ordered phase and NSL phase the transition is of the 
first order. As the value of /3 is further increased, the intralayer correlation becomes weaker 
and the second order transition to the ILD phase takes place at /? = 4. This is verified 
analytically from the stabilty analysis of Eqs. ( p.44| , pl45| , p.46|) around the ILD solution. Of 
course, this phase transition between the NSL phase and ILD phase is an artifact of the 
mean field approximation and it should be interpreted as a crossover. For very small or 
large value of a, the direct first order transition from the Neel or collinear ordered phase to 



the ILD phase takes place as in the unfrustrated case. pO| 

It should be remarked that the spin-fluid phase with considerable intralayer correlation 
is stabilized around a ~ 0.5 where the effect of the frustration is most pronounced. In this 
sense, our NSL state is the first promising example of the non-trivial frustration induced 
quantum spin fluid phase in two dimensions. 

Figures Ka) and ||(b) show the amplitude of the sublattice magnetization and Nq 
along the line with fixed a. In general, the sublattice magnetization is enhanced for small (3 
and turns to decrease for larger f3. Namely, the interlayer coupling strengthens the ordering 



as far as it is small. The same is true also in the unfrustrated case. 20,21 This feature 



can be already observed within the linear spin wave analysis. Dotsenko pO[ has also 



obtained the similar result using the mapping onto the nonlinear cr-model. 
The energy gap in the NLS phase AE is given by 



AE = E+{k = 0) = ^(8A + + - 5). (4.54) 

Figure |^ shows the variation of AE for various values of a. It grows from starting from the 
Neel-NSL boundary and saturates at the NSL-ILD boundary. At the coUinear-NSL boundary 
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the energy gap in the NSL phase remains finite reflecting the fact that this transition is of the 
first order. However, we may expect that this is also an artifact of the present approximation. 
Taking into account that AE is proportional to the inverse of the spin-spin correlation length 
C, for ^ >> 1, the antiferromagnetic short range correlation is highly enhanced in the NSL 
phase near the phase boundary. 

V. SUMMARY AND DISCUSSION 

The spin-1/2 bilayer Ji — J2 model is studied by means of the modified spin wave ap- 
proximation and the ground state phase diagram is obtained. For small interlayer coupling, 
Neel or coUinear type long range order exists for any value of a. However, with the increase 
of the interlayer antiferromagnetic coupling f3, the correlated quantum spin fiuid phase ap- 
pears around a ~ 0.5. The width of the spin fiuid phase becomes wider as f3 increases. On 
the other hand, the correlation length decreases with the increase of (3 and we only find the 
interlayer dimer phase for j3 > 4. 

It has been widely expected that the frustation effect enhance the quantum fiucutation 
and leads to the quantum spin liquid phase in two dimensional system. Although the single 
layer Ji — J2 model is one of such candidates, the conclusion is rather sensitive to the 
approximations used, the method of numerical calculations and data analysis. On the other 
hand, we may expect the presence of the frustration induced highly correlated quantum spin 
liquid over a wide range of parameters for the bilayer Ji — J2 model. 

We have also found that the Neel state remains stable for the value of a slightly larger 
than 0.5. This may be explained as follows: In the collinear phase, the classical ground state 
is continuously degenerate and therefore the quantum fiuctuation is more pronounced than 
the Neel phase. Therefore the Neel state is stabilized rather than the collinear phase even 
for a;>0.5. 

Using the modified spin wave approximation, Nishimori and Saika ^ obtained the result 
that the ground state energy jumps at the transition from the Neel phase to the collinear 
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phase in the spm-1/2 single layer Ji — J2 model, while in our calculation this transition is 
a usual first order transition even for (3 = 0. This is due to the fact that Nishimori and 
Saika expanded the ground state energy with respect to l/S* and trunciated at the second 
order. Actually, such estimation is known to give the results better than the mean field 



type estimation in some cases. pil| However, here we employ the naive mean field ground 



state energy without expansion for the consistency within the present approximation. In 
any case, such details of the transition are beyond the scope of our approximation. 

Althouth we expect that our approach captures the essential features of the ground state 
of the present model, our approximation is far from quantitative. Even in the unfrustrated 
case {a = 0), modified spin wave theory predicts rather large Neel-ILS critical value of 
/3 = /3c ~ 4.25 compared to the more reliable estimation /^^ ~ 2.5 by the dimer expansion 



!T|] and the quantum Monte Carlo simulation. Unfortunately, the quantum Monte Carlo 
simulation is not expected to be powerful enough for the frustrated model due to the negative 
sign problem. The dimer expansion method may be promising but higher order calculation 
requires excessive computational time and memory. Further investigation is thus required 
for the quantitative understanding of the present model. 
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work is supported by the Grant-in-Aid for Scientific Research from the Ministry of Educa- 
tion, Science and Culture. The numerical calculation is performed by the HITAC S820/15 
at the Information Processing Center of Saitama University. 
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FIGURES 

FIG. 1. The ground state phase diagram of the bilayer Ji — J2 model. The sohd and broken 

lines are the lines of the first and second order transition, respectively. 

FIG. 2. The /3-dcpcndence of the sublattice magnetization in the (a) Neel and (b) collinear 
ordered phases. The values of a are indicated in the figure. The values on the phase boundary are 
represented by the open circles. 

FIG. 3. The /3-dependence of the excitation gap in the disordered phase. The values of a are 
indicated in the figure. The values on the phase boundary are represented by the open circles. 
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